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WAVE SURFACES DUE TO IMPACT ON ANISOTROPIC FIBER 
COMPOSITE PLATES 


by Francis C. Moon* 

Lewis Research Center 

SUMMARY 

The use of advanced fiber composite materials for fabricating fan or compressor 
blades depends on the ability of the materials to sustain limited damage under impact 
forces. For foreign objects whose masses are comparable to that of the blade, such as 
large birds, one can expect impact times of the order of the lower blade vibration per- 
iods. For these cases the method of vibrations should prove adequate to the task of pre- 
dicting the impact stresses and possible failure modes. For smaller objects such as peb- 
bles or hailstones, the impact times are of the order of microseconds. Thus the energy 
transmission to the blade is completed before there are any or many reflections from the 
boundaries. It is for this class of problems that a stress wave analysis is more useful 
than vibratory methods of analysis. The stress waves induced in anisotropic plates by 
transverse, short-duration impact forces are examined in this report. The anisotropy 
is related to the layup angles of the fibers, which lie in the plane of the plate. Using a 
modification of Mindlin's approximate theory of plates, it is shown that both extensional 
and bending waves are generated by transverse impact. The magnitudes of the wave 
velocities in different directions are calculated for graphite fiber -epoxy matrix plates as 
well as boron -aluminum and glass -epoxy systems for various layup angles. Finally, the 
shapes of the wave fronts or wave surfaces due to point impact are also presented for 
the cases mentioned. 


INTRODUCTION 

The successful application of advanced fiber composite materials to jet engine fan or 
compressor blades will depend in part on the ability of these materials to withstand the 
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forces of impact due to foreign objects. Such impact can be the result of the ingestion 
of stones, nuts and bolts, hailstones, or birds into a jet engine. The relative velocity 
of the impacting body to the blade can be in the order of 450 meters per second (1500 ft/ 
sec). The ingestion of objects of sizeable mass (e.g. , birds) might involve the dynamics 
of the entire blade. The high speed impact of small objects will result in small impact 
times (<50 psec), and the initial transmission of impact energy into a local region of the 
blade. This initial energy will propagate into the rest of the blade in the form of stress 
waves . Although such high speed impact will involve local cratering or even complete 
penetration, long range damage away from the impact area can result from the reflection 
of stress waves (spalling) and focusing due to changes in geometry. 

It is also observed that materials under high rates of strain exhibit an increased ten- 
sile strength and decreased ductility. Such evidence tends to validate the use of elastic 
wave analysis for the determination of the prefracture stresses induced by the impact 
forces. However, even if plastic waves do predominate, elastic precursor waves will 
bound the stressed impact zone. 

In this report, calculations of velocity and wave surfaces in anisotropic composite 
plates due to transverse impact forces are presented. These wave surfaces, for a given 
time after impact, bound the stressed region surrounding the impact point. 


REVIEW OF BASIC EQUATIONS OF ANISOTROPIC ELASTIC PLATES 


The composite plates under consideration are imagined to comprise a number of uni- 
directional plys (fig. 1). An equal number of plys lie at angles ±<p from the symmetry 
axis in such a manner that bending -extensional coupling does not result. We also assume 
that the number of plys across the plate thickness is reasonably large, so that average 
properties across the plys can be used. This approximation will be valid for wavelengths 
greater than the ply thickness and certainly valid for wavelengths greater than the plate 
thickness. 

Thus, in place of the n-ply plate, the wave propagation in an equivalent anisotropic 
plate is being studied. The equivalent elastic constants are obtained from a static anal- 
ysis of the n-ply composite plate. 

The equations of motion for a linear anisotropic elastic body are (ref. 1) 



(1) 


(The double summation convention is assumed, the dots indicate time differentiation, and 
the notation cp . stands for dcp/dx. where x. is a Cartesian coordinate.) The vector & 
is the displacement, and body forces are assumed to be absent. In general, the stress 
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tensor is related to the strains e^ through the equations (ref. 1) 


*1] “ C ijkZ e kZ 


( 2 ) 


There are only 21 independent elastic constants in general. For orthotropic sym- 

metry, which pertains to the composite plates tinder consideration, the stress-strain 
equations take the following matrix form: 
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of course related to the 


are related to the displacements in the body by 

e i] = r (u i,j + u i,i ) 


(3) 


(4) 


Combining equations (1), (2), and (4) results in the following partial differential equa- 
tions: 
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ijkZ U k , £ j 
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Wave propagation in anisotropic media has been studied for a long time (refs. 2 and 3); 
however, very few problems have been solved in which boundaries are present. 

The approximate theory of anisotropic plates to be used in this study is due to Mindlin 
and coworkers (refs. 4 and 5). In their theory the three-dimensional displacement is ex- 
panded in Legendre polynomials in the thickness direction. 
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where 17 = Xg/b and b is half the plate thickness and where 

P 0 fa) = 1 

PjM - V 

p M = - 3 t? 2 - 1 
z 2 

The { P n (?7) } are orthogonal: 



The functions jjf 1 ^ have a physical significance (see ref. 5, pp. 563-564): u!?,u2 repre- 

0 1 « ' 

sent in-plane or extensional deformation; u„ represents the transverse displacement of 
11 “ 

the plate; Uj and are measures of the bending strains or b b ^ where i//j is 

the slope of the plate midsurface due to bending about the 3 -axis; and u^ is a measure 
of the thickness stretching. 

To obtain the approximate equations of motion, a variational method is used 
(ref. 4). Instead of solving equation (5) directly, the equations are integrated across the 
thickness: 



- pii.)5u.b dri dA = 0 
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where 5m and t^. are calculated using the series representation of the displacement 
(eq. (6)). This leads to an infinite set of equations each involving higher modes of vibra- 
tion of the plate: 
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where a = 1, 3. 
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If impact forces are present on the upper surface of the plate, then the following boundary 


conditions are used to evaluate the second terms, 



in equation (8) . 
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t 21 (i? = ±1) = 0 
t 23 (rj = ±1) = 0 
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This scheme has been carried out for n = 0, 1, 2 for orthogonal symmetry. The 
equations of motion for ulj* and u 3 and u 2 , u|, and u 3 are shown in equations (10) 
and (11), respectively: 
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The first two equations describe the in-plane or extensional motion. Note that this 
motion depends on the thickness stretching Ug. The equation of motion for u*, which is 
obtained from equation (8) , is 
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The following approximation is made. We drop the higher order displacements in the 
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displacement expansion, Uj , u^ , etc. , in equation (6) and in equation (12). Next the 
terms containing second derivatives are dropped in equation (12), keeping only the low- 

frequency terms. This procedure leaves an explicit equation for ui. A similar pro- 
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cedure is used 


sed on the equation for u), ' (which is not presented in full here). The sim- 
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plified equations for u 2 , u), ' are 
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Using these equations, the terms u* and can be eliminated from equations (10) 
and (11). The resulting set of equations form the basis of our wave analysis: 
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It should be noted that in the procedure used by Mindlin (ref. 4) the coefficients C 44 
and Cgg in equations (10) and (11) were replaced by k g C 44 and kjCgg, respectively. 
The correction constants and kg were adjusted in order to match the thickness 
shear vibration mode (ref. 4). These terms will not enter the calculations presented 
herein. 


WAVE PROPAGATION 

Solutions to the propagation of plane waves in anisotropic plates are sought in the 
form 


& = ,&&• vt) (16) 

where ji is a prescribed unit vector. If such a solution exists, & can only change in the 
direction jj,» that is, 
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or at a given time the displacement is constant on a plane normal to the vector ji. When 
a solution of the form (16) exists, v is called the wave or phase speed. 

Consider first the extensional motion which is governed by equations (14). Assume a 
solution in the form 


uj = Ujf(n, • x - vt) 
u 3 = u 3 f( a • & - vt) 


(17) 


Substituting these expressions into equations (14) reveals that Uj, Ug, and v must sat- 
isfy the following linear algebraic equations for a given n: 
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Thus v is a root of the equation 


A(a, v 2 ) =det^A.. - pv 2 6 ij ) = 0 


(19) 


where is the Kronecker delta = ®21 ~ ® and ^11 = ®22 = Physically 

possible elastic constants will guarantee that A„ is positive definite. This guar- 

2 2 

antees two positive real roots and Vg for a given wave normal ji. 

The ratio U^/Ug will be determined by substituting each root v^ into the equa- 
tions (19). Since A„ = A^, the displacement vectors corresponding to the roots v^ 
and Vg will be orthogonal to each other. If the displacement direction, determined by 
Uj/Ug, is parallel to ji , the wave is called longitudinal; if the displacement correspond- 
ing to Uj/Ug is normal to &, the wave is called transverse. For isotropic materials it 
is known that the wave motion is longitudinal for the larger root, and transverse for the 
smaller root. For anisotropic materials, the velocity v depends on ji , and the motion 
is neither longitudinal nor transverse except for certain symmetry directions. 

Consider next the bending equations (15). One can show that the only plane wave so- 
lutions of the form (16) that satisfy equations (15) are harmonic functions, that is, 
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The product w = kv is called the frequency; k is called the wave number or inverse 
wavelength. For bending motion, the phase velocity v depends on the frequency u> as 
well as the wave normal ji. Mindlin (ref. 4) has examined the dependence of v on w 
for various material anisotropies. 

Thus the behavior of the bending motion at the wave fronts cannot be determined in 
the same manner as was the extensional motion. Instead of finding a valid solution for 
the whole impact disturbed area of the plate, we consider the motion at the wave front 
only. Across this front, one imagines that certain quantities have discontinuities. The 
displacement and the stress are assumed to be continuous across the wave front but allow 
discontinuities in the second derivatives of U. Such waves are called acceleration waves 
(ref. 6, chapt. 5). 

Let [^/] denote the jump in the function \p(x^,x^) across the wave front. Then by as- 
sumption we have (where i, j = 1 or 3) 



The second derivatives are assumed to exist on both sides of the wave front; thus, we can 
write the equations of motion for bending (eq. (15)) on both sides of the wave front and 
subtract one from the other , from which results 
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The jump in acceleration, however, is not independent of the jump in the strain 
gradient. It can be shown that for a plane wave front with unit normal the following 
relations hold: 
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The quantity v is called the wave front speed in the normal direction. This relation can 
then be used in the preceding equations to obtain linear algebraic relati between the 
discontinuities in acceleration across the front: 
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and are exactly the same constants that occur in equation (18). 

Thus the wave fronts associated with a jump in the bending accelerations 
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2 1/2 

and 9 u Q /3t travel at the same speeds as the wave front associated with the extensional 

«/ 2 0/2 
motion. There is another wave front corresponding to a jump in the quantity 8 ulj /8t . 

For the case of a composite with symmetric ply orientation about the midplane, 

C 66 = C 44 

0/2 

The bending wave associated with the jump 8Ug /9t is thus isotropic. 

If both extensional and bending motions are generated simultaneously by impact, the 
two extensional and two bending wave fronts will travel with the same wave speeds . 

The analysis presented here is not unique. The same results can be obtained if one 
considers the equations of motion (14) and (15) from the method of characteristics 
(ref. 6). 

The velocity surfaces v = v(jj) have been computed for various fiber composites and 
ply configurations. These results are discussed in a later section. 

WAVE SURFACES 

In the preceding section we outlined how plane waves would travel in an anisotropic 
plate. The phase velocity of two of the modes was found to depend on the orientation of 
the wave normal to the symmetry axes of the plate. This angle we called cp. Suppose, 
then, that a plate receives a transverse impact at the origin of a coordinate system 
(r, 6) . The disturbance can be thought of as a superposition of plane waves. To an ob- 
server at position (fq, 9q), the first signal to arrive may not be that corresponding to the 
wave normal cp = 6 q. If t is the arrival time, the first plane wave n(<p) to arrive at 
the point £ must satisfy (see fig. 1) 


£ • n((p) = v(<p)t 

For a given time (say t = 1) the wave surface is defined as the locus of points £, which 
satisfy (unpublished notes by Yih-Hsing Pao) 
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The vector g, is called the slowness vector, and the surface l/v(<p) is called the slow- 
ness surface. (A good discussion of the properties of velocity, slowness, and wave sur- 
faces may be found in ref. 2.) 

Instead of finding the first arrival wave &(<p) for a given r , we determine £ for a 
given plane wave jj, such that equation (24) is satisfied, with £ fixed. Then, the equa- 
tion g, • £,= 1 represents a line in the slowness plane and £ the normal to that line. 
However, not all g, are admissible; g, has to be on the slowness surface. Thus, the 
line g, ■ £,= 1 is tangent to the slowness surface and £ is the normal vector to that 
surface . Suppose the slowness surface is given by the equation 

-g(a) = 0 (25) 

Then 
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Substituting this expression for r into equation (24) yields 
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In our case v and hence 1/v are roots of a quadratic equation (19), that is, 
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For each root v there is a wave surface. It can be shown (ref. 3) that the outer 
surface, which is associated with the fastest velocity, is strictly convex. However, the 
slower velocity surface can result in a wave surface with cusp points. 

The locus of Jiicp) has been computed for various fiber composite systems and fiber 
layup angles. The results are shown in figure 2 and are discussed in the following sec- 
tion. 


DISCUSSION OF NUMERICAL RESULTS 

Velocity, slowness, and wave surfaces were calculated for various anisotropies 
corresponding to various fiber composite plates using a digital computer. The three 
fiber-matrix systems examined were graphite -epoxy, boron -aluminum, and glass-epoxy. 
The equivalent elastic constants for these fiber-matrix systems at various layup angles 
(fig. 1) were obtained by Chamis (ref. 7). These constants, which are listed in tables I 
to m, are based on a statistical analysis of an eight-ply plate using the known properties 
of each fiber -matrix ply. 

The velocity and wave surfaces for a boron -aluminum composite are shown in fig- 
ure 2. The ratio of moduli Cjq/Cgg = 1.2. This results in different wave speeds in the 
two principal directions. However, the shear velocity is almost isotropic. The wave 
surfaces (fig. 2(b)), while showing the effects of anisotropy, exhibit no peculiarities. 

The graphite -epoxy system contrasts with the boron-aluminum system because of its 
high stiffness ratio; Cqj/Cgg = 24 (zero layup angle). The velocity surfaces for layup 
angles of ±0°, ±15°, ±30°, and ±45° are shown in figure 3. It is interesting to note that, 
as the fiber orientation approaches ±45°, the anisotropy in the larger wave velocity 
(quasi-longitudinal wave) diminishes, but that of the smaller root (quasi-shear wave) in- 
creases. 

The resulting wave surfaces for graphite- epoxy are shown in figure 4. (The slow- 
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TABLE I. - STRESS-STRAIN COEFFICIENTS FOR 55 PERCENT BORON 
FIBER -ALUMINUM MATRIX COMPOSITE 

[Material density, 2. 65 g/cm 3; all constants to be multiplied by 10 6 psi; data obtained from ref. 7.] 




0° Layup 





15° 

Layup 




42.80 

11.44 

11.44 

0 

0 

0 



42.14 

11.67 

11.54 

0 

0 

0 




34.47 

14.92 

0 

0 

0 




34.47 

14.68 

0 

0 

0 





34.47 

0 

0 

0 





34.92 

0 

0 

0 






9.775 

0 

0 






9.775 

0 

0 







13.18 

0 







13.29 

0 








9.775 








9.775 





±30° 

Layup 





±45° 

Layup 




40.40 

12.31 

11.75 

0 

0 

0 



38.22 

13.18 

11.85 

0 

0 

0 




34.47 

14.05 

0 

0 

0 




34.47 

13.18 

0 

0 

0 





36.24 

0 

0 

0 





38.22 

0 

0 

0 






9.775 

0 

0 






9.775 

0 

0 







13.49 

0 







13.60 

0 








9.775 








9.775 



TABLE II. - STRESS-STRAIN COEFFICIENTS FOR 55 PERCENT GRAPHITE 


FIBER-EPOXY MATRIX COMPOSITE 

[Material density, 1. 44 g/ cm 3 ; all constants to be multiplied by 10 6 psi; data obtained from ref. 7.] 


0° Layup 


27.95 0.3957 

1.170 


0.3957 0 0 

0.4601 0 0 

1.170 0 0 

0.3552 0 

0.7197 


±30° Layup 


16.48 0.4118 5.167 0 

1.170 0.4400 0 

3.093 0 

0.3552 


0 

0 

0 

0 

5.491 



±15° Layup 

0 

0 

0 

0 

0 

0.3552 



24.56 0.4000 1.986 0 0 0 

1.170 0.4558 000 

1.374 0 0 0 

0.3552 0 0 

2.310 0 

0.3552 



±45° Layup 

0 



8.197 0.4279 6.758 000 


0 



1.170 0.4279 000 


0 



8. 197 0 0 0 


0 



0.3552 0 0 


0 



7. 082 0 


0.3552 



0.3552 
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TABLE m. - STRESS-STRAIN COEFFICIENTS FOR 55 PERCENT GLASS FIBER-EPOXY 

MATRIX COMPOSITE 

Q fi 

[Material density, 1. 92 g/cm ; all constants to be multiplied by 10° psi; data obtained from ref. 7.] 




0° Layup 





±15° Layup 





7.500 

0.9097 

0.9097 0 

0 

X) 



6. 890 

0.9321 

1.178 0 

0 

0 




2.395 

1.244 0 

0 

0 




2.395 

1.222 0 

0 

0 





2.395 0 

0 

0 





2. 468 0 

0 

0 





0. 5753 

0 

0 





0. 5753 

0 

0 






0.9457 

0 






1.214 

0 







0. 57 53 







0. 5753 




±30° Layup 





±45° Layup 





5. 419 

0.9933 

1.715 0 

0 

0 



3. 874 

1.077 

1.983 0 

0 

0 




2. 395 

1. 161 0 

0 

0 




2.395 

1.077 0 

0 

0 





2. 866 0 

0 

0 





3. 874 0 

0 

0 





0. 5753 

0 

0 





0. 5753 

0 

0 






1.751 

0 






2.019 

0 







0. 5753 







0. 5753 



ness surface is shown in fig. 5.) In contrast to the boron -aluminum system, the quasi- 
shear surfaces show peculiar cusps and nonconvexity. This behavior is also characteris- 
tic of crystal systems such as zinc. Unlike the natural crystals, we can change the wave 
properties, without changing the material constituents, by varying the fiber layup angle. 

It becomes clear that, as the anisotropy in the quasi-longitudinal wave is reduced (i.e. , 
cp — 45°), the cusped behavior of the quasi-shear waves increases. This is due to the 
previously mentioned increase in shear wave anisotropy as cp — 45° (fig. 3). 

Another peculiar property of wave propagation in this composite system can be noted 
by examination of the ±45° fiber layup case (fig. 4(d)). On the outer wave surface, the 
angle of the wave normal of the first arrival plane wave is listed. One can see that the 
distribution of plane wave normals is heavily concentrated at positions on the wave sur- 
face close to the fiber directions. This might imply a focusing of waves along the fiber 
directions. For the other fiber orientations, the distribution of wave normals is also 
concentrated at those points on the wave surface close to the fiber directions but not as 
densely as in the ±45° layup case. The implications of this wave focusing along the fiber 
direction will not be made completely clear until the stress and displacement fields are 
found. 
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Similar results for the glass fiber-epoxy composite system are presented in figures 
6 and 7. The ratio of stiffnesses for this case is = 3.1 (zero layup angle). The 

wave surfaces for this system show features similar to the graphite-epoxy case. Note, 
however, that for a layup angle of ±15°, the quasi-shear wave velocity is almost isotropic 
(fig. 6). This results in a wave surface (fig. 7(c)) -with no cusped behavior. Although not 
as marked as the graphite -epoxy case, this system also exhibits a wave normal focusing 
along the fiber directions. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, March 5, 1971, 

129-03. 
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Figure 1. - Diagram notation and sign convention for 
plate element. 
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(a) Velocity surfaces. 


Distance from 
impact point, 
mm 



(b)Wave surfaces at 1 microsecond after impact. 


Figure 2. - Velocity and wave surfaces for 55 percent boron 
fiber - aluminum matrix plates. Fiber layup angle, 0°. 



Velocity, 

mm/Msec 


Figure 3. - Velocity surfaces for 55 percent 
graphite fiber-epoxy matrix plates at various 
layup angles. 
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(a) Fiber layup angle, 0°. (b) Fiber layup angle, ±15°. 


Figure 4. - Wave surfaces at 1 microsecond after impact for 55 percent graphite fiber - epoxy matrix plates. 
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(c) Fiber layup angle, +30°. 
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Figure 4. •'Concluded. 


(d) Fiber layup angle, ±45°. 




Figure 5. - Normalized slowness surfaces for 55 percent graphite 
fiber - epoxy matrix plates at various fiber layup angles. 


22 


ii ■lamiin it 


111 lllllll III II 


Distance from 
impact point, 
mm 


Distance from 
impact point, 
mm 


(a) Fiber layup angle, 0°. (b) Fiber layup angle, ±15°. 

Figure 7. - Wave surfaces at 1. 15 microseconds after impact for 55 percent glass fiber - epoxy matrix plates. 
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(c) Fiber layup angle, ±30°. (d) Fiber layup angle, ±45°. 

Figure 7. -Concluded. 
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